Structural bioinformatics analysis of SARS-CoV-2 variants reveals higher hACE2 receptor binding affinity for Omicron B.1.1.529 spike RBD compared to wild type reference

To date, more than 263 million people have been infected with SARS-CoV-2 during the COVID-19 pandemic. In many countries, the global spread occurred in multiple pandemic waves characterized by the emergence of new SARS-CoV-2 variants. Here we report a sequence and structural-bioinformatics analysis to estimate the effects of amino acid substitutions on the affinity of the SARS-CoV-2 spike receptor binding domain (RBD) to the human receptor hACE2. This is done through qualitative electrostatics and hydrophobicity analysis as well as molecular dynamics simulations used to develop a high-precision empirical scoring function (ESF) closely related to the linear interaction energy method and calibrated on a large set of experimental binding energies. For the latest variant of concern (VOC), B.1.1.529 Omicron, our Halo difference point cloud studies reveal the largest impact on the RBD binding interface compared to all other VOC. Moreover, according to our ESF model, Omicron achieves a much higher ACE2 binding affinity than the wild type and, in particular, the highest among all VOCs except Alpha and thus requires special attention and monitoring.


Results
In order to enable a comprehensive structural analysis of emerging partial or complete SARS-CoV-2 genomes, we follow a three-phase approach. In the first phase, we analyze new genomes (e.g. from GISAID or sequences directly provided by associated laboratories) for sequences with amino acid exchanges in regions of SARS-CoV-2 protein structures that are potentially relevant for drug binding or cell uptake.
In phase two, we look for sequences where one or more mutations are found that have not yet been investigated in our structural analysis pipeline. In such a case, the sequence is submitted to a structure modeling workflow, whose results can be used to predict an influence caused by the respective amino acid exchanges on active site regions or cavities that are potentially relevant for drug development. In the case of the spike protein, on whose RBD we focus herein, relevant models are analyzed using our Catalophore Halo technology (defined and illustrated below).
In phase three, we compare any modified RBD's Halo to the wild type Halo. If such a Halo comparison shows a substantial change, we expose the new variant to an ESF-based molecular dynamics modeling pipeline in order to predict the corresponding change in binding affinity.
SARS-CoV-2 sequence analysis of spike-RBD diversity. By December 6th 2021, 5,799,116 SARS-CoV-2 genome sequences and 5,649,261 spike-protein sequences were available at GISAID 21 . For the sequence analysis of the SARS-CoV-2 RBD diversity we extracted 5,173,253 RBD sequences by filtering the spike protein sequences for the RBD flanking residues "TSNF" (position 315-318) and "NFNG" (position 542-545) as well as for correct size of 223 amino acids. Within these sequences we identified a total number of 7,700,325 amino acid exchanges, which amounts to 0.67% of all residues. Compared to August 2020, when 185 mutations were identified in 73,042 spike RBD sequences 13 , this occurrence per sequence of mutations shows a 588-fold increase, highlighting the progressive evolution of this virus 30,31 . The accumulation of mutations was significantly higher Visualization of amino acid exchanges using difference Halos. The influence of respective amino acid exchanges in the spike RBDs of the VOCs Alpha, Beta, Delta and Omicron regarding crucial properties can be visualized by Catalophore Halos. Halos are equidistant point clouds with values representing 16 physicochemical property fields such as electrostatics, hydrophobicity, flexibility, potential energies, hydrogen-bonding acceptor/donor potentials, dissolvability, aromaticity, and others in a discretized manner, surrounding the molecular surface of the molecule that induces these fields 39 . Originally, they were developed for the calculation of fields inside buried protein cavities 40 . In the case of spike-RBD-hACE2 complexes, the Halo point clouds are   (Fig. 2b,c). Interestingly, the Halos of Beta and Gamma appear to be  www.nature.com/scientificreports/ more different from the wild type regarding electrostatics than Delta (Fig. 2a), whereas Alpha shows almost no difference regarding electrostatics. Regarding hydrophobicity, Alpha reveals a local increase of hydrophobicity (blue patch) in the region around its N501Y mutation compared to the wild type. Due to the recurrence of the same blue patch in the difference Halos of Beta, Gamma and Omicron compared to the wild type, it may be linked to alterations in biomolecular interactions and should be investigated in further detail. By visual comparisons of electrostatics and hydrophobicity difference Halos, Omicron seems to be akin to a combination of Beta and Delta with additional changes in the region where mutation N440K is located (Fig. 2b). The large Halo difference between Omicron and Delta is explained by mutations G339D, S371L, S373P, S375F, K417N, N440K, G446S, S477N, E484A, Q493R, G496S, Q498R, N501Y, Y505H that are only present in Omicron and the absence of L452R while only sharing T478K. At the binding interface, the new arginine replacing glutamine at position 493 in Omicron yields a new strong hydrogen bond with hACE2 amino acid E35 which is clearly reflected by the dark blue area (increase in positive charge) of the Omicron-WT difference Halo as illustrated by Fig. 2c. These depictions allow a rapid qualitative overview of the extent of physico-chemical changes at the RBD-hACE2 binding region. As a consequence, they allow a quick decision on which variants to prioritize for further resource-intensive structural and dynamics analysis, as for instance calculating binding affinities to hACE2 through the ESF we developed for this purpose.
Development and optimization of a binding affinity model. We employed a two-parameter ESF closely related to the LIE method 41,42 where the parameters are associated with van der Waals (VDW) and electrostatic interaction differences between the bound and unbound state. While different sets of physical assumptions lead to concrete values for these parameters in the original LIE method, they can also be used, based on a set of experimental data, to gauge the contribution of the corresponding terms to changes in ΔG [43][44][45] . Further details on the ESF setup as well as on various model assumptions and tests can be found in the methodology section below and in the supporting information.
In essence, our approach for a two-parameter model turned out to be a reasonable starting point as well as the most sensible model regarding the number of parameters. We confirmed this by a fit (for each case) to a suitable set of experimental data (as detailed below) as well as a round of leave-one-out cross-validation calculations that predict a set of ΔG. The correlation of this prediction with the experimental set is our concluding point of reference for model usefulness. As the analysis below shows, the resulting weight of the electrostatic contribution is very small and, within errors, consistent with zero. This may seem surprising, given the physical origin and motivation of this term, but straight-forwardly results from the fact that Coulombic energy difference values exhibit statistical errors of the same magnitude as the values themselves. The fitting procedure translates this weak signal-to-noise ratio into a very small weight for the electrostatic term.
The most obvious explanation for the large statistical width of electrostatic contributions is connected to the huge binding interface of the hACE2 ligand comprising around 35 amino acids (ca. 500 atoms), compared to simple protein-ligand interfaces with much smaller ligands and therefore much less noise. In addition, the Coulomb potential acts over longer ranges (compared to the Lennard-Jones potential which quantifies energy contributions due to short-range VDW forces) and its evaluation involves substantially more interacting atom pairs including many highly flexible surface side chains. The more rapidly declining VDW interactions appear to be more stable and significantly less sensitive to small displacements in atom positions 13 and interface size. By design, the intermolecular Lennard-Jones potential is mostly affected by a narrow region of the protein-protein interface involving far fewer atoms than in case of long-range electrostatics. As a consequence, VDW contributions emerge as a clear signal and furthermore correlate well with experimental data. This observation indicates a strong dependence of the binding affinity on the steric fit of the two involved proteins. Thus, we are able to pioneer the application of an LIE related empirical scoring function to absolute protein-protein binding affinity prediction. Other published studies of LIE applied to protein-protein binding do exist but were limited to local subdomain analysis as for instance of regions around single mutations 46,47 .
Despite this good performance of a model with (effectively) a single parameter, we refrain from dropping the electrostatic term completely for the sake of generality. It is unclear, from our study alone, whether or not the term would be beneficial and/or significant in other cases.
According to the last prediction column in Table S1 (visualized in Fig. 3a), the highest prediction accuracy and lowest error was achieved with 50 replicates of 200 ps molecular dynamics (MD) simulation trajectories, where, for model fitting and cross-validation respectively, R 2 amounts to 0.77 and 0.74 with an average error of less than 1.9 and 2.0 kJ/mol.
A corresponding scatter plot of fitted and cross-validated vs. experimental binding energies ΔΔG (relative to the wild type WT) is depicted in Fig. 3c. Fortunately, the experimental ranking of our training set variants is very well reflected by these model predictions. In addition, the distribution of experimental ΔG values in terms of mean and standard deviation (− 55.76 ± 4.7 kJ/mol) agrees extraordinarily well with those of the fitting and cross-validation procedure, both amounting to − 55.8 ± 3.9 kJ/mol. Furthermore, overfitting can safely be considered as minimal since corresponding points (variants) of fitting and cross-validation plots as well as both regression lines are lying very nearby. The results in Fig. 3a have been obtained with the optimal fitted weights, w vdw = 0.76472 and w elec = 0.02355, associated with model A.1. Having applied this predictive model to an unpublished internal set of 21 hACE2 variants in complex with the wild type for which EC50 values were available, we obtained a remarkable prediction accuracy of R 2 = 0.61. Prior to that, EC50 values had been converted to binding free energies in analogy to K D values.
Considering the number of replicates needed for a sufficient convergence of the energies and consequently of affinity-based variant rankings, the plot in Fig. 3b indicates moderate stability already from 30 replicates on, although MAE and R 2 are still slightly worse compared to 50 replicates (Tables S2-S4). We can safely assume a  Table S1). However, with regard to the required run time (currently 1.1 days on a Amazon AWS C6i.16xlarge instance in case of 200 ps trajectories), 50 replicates seem to be a reasonable choice if we put more emphasis on accuracy. For this convergence analysis we determined the fraction of variant pairs associated with a swap in their accumulative energy order while increasing the number of replicates. All frequencies in percent were obtained via division by the total number of variant pairs n = 43 2 = 1849. In more urgent cases, 30 replicates would take around 15 h to run at a slightly lower accuracy. Another major quality aspect of binding affinity models is related to the predicted ranking of VOCs, as we strive for high consensus between predicted and experimental top-N variants. Figure 3d illustrates the fraction of variants among the top N predicted ones that are also included in the top N experimental variants of the training set. We compared the performance of an empirical model (a) trained with 50 replicates and applied to predictions on the basis of 50 replicates as well, (b) same as (a) but using 30 replicates for both steps, and (c) a mix of (a) and (b), namely model training with 50 but predictions on 30 replicates. Finally, (d) a non-empirical prediction model published by Singh et al. 13 in 2020 and based on umbrella sampling along with weighted-histogram analysis was applied to the same training set of 43 variants. As expected, the highest consensus, especially in the range of up to top 5, was obtained by models trained with 50 replicates, namely cases (a) and (c) corresponding to the green and red plots (where green is sometimes hidden behind the red line). Regarding predictions of training set items, up to 80-100% of the top 3-5 are in agreement with experimental top 3-5 candidates according to the ESF approaches (a) and (b) and significantly less consensus along the entire range of N variants using the umbrella sampling method, in particular less than 50% among the top 5 variants.
Affinity model application to VOCs. Using our final predictive model calibrated with the optimal weights, we estimated binding affinities for VOCs that had emerged during the past year: Alpha, Beta, Gamma, Delta, and most recently Omicron. However, we have to bear in mind that our structural variant models only comprise the spike RBD rather than the entire spike trimer. Table 1 compares our calculated relative binding energies ΔΔG  According to the results in Table 1 and Fig. 3c, the predicted order of dissociation constants for WT, Alpha, Beta and Gamma exactly reflects experimental findings 32,48 and yields a coefficient of determination R 2 of 0.88 based on ΔΔG values. In particular, with a 14-fold decrease compared to the wild type dissociation constant, the most outstanding binding affinity among all VOCs was predicted for Alpha. This is also in line with lab data 32,48 revealing around tenfold binding affinity increase and reports 49 stating, for the UK between October and November 2020, a 70-80% rise in transmissibility with respect to the wild type reference.
For the new SARS-CoV-2 Omicron spike variant, our studies reveal a similarly remarkable 7.7-fold increase of binding affinity compared to the wild type. In terms of binding strength it therefore ranges between Alpha and Gamma. A significantly lower binding affinity, somewhere around Beta but lower than the wild type, was predicted for Delta and Delta plus. Moreover, Delta plus (having K417N in addition to the amino acid exchanges found in Delta) was predicted to have a lower affinity than Delta, which is supported by experimental results revealing weaker binding due to the presence of K417N 48,50 . In addition, the authors allude to an increase of immune evasion caused by this mutation. It should be noted that Omicron as well includes this mutation possibly explaining its lower binding affinity compared to Alpha and thereby indicating a strong tendency to immune evasion for Omicron. The relative order of WT, Delta and Omicron is also in line with Docking results recently published by Kumar et al. 51 , although energy magnitudes and therefore differences given in that article translate to a K D ratio of 10 million for Omicron compared to the wild type.

Discussion
Following the trade-off model of virulence 52 , the SARS-CoV-2 virus, like other viruses 53 , is constantly evolving 54 by modulating the rate of infectious transmission, higher virulence, and higher virus production in order to improve viral fitness. In this context, viral immune evasion is an evolutionary strategy to allow for the coexistence of viruses and their hosts 55 as shown with artificial polymutant SARS-CoV-2 spike protein pseudotypes that resisted antibody neutralization to a similar degree as circulating VOCs 56 .
Typically, viruses adapt to a specific host. However, when hosts fluctuate in time or space, generalist viruses may evolve as well 57 . It is hypothesized that during this process, intermediate virulence maximizes pathogen fitness as a result of a trade-off between virulence and transmission 58 . Although it is often speculated that virulence of viruses will finally decrease over time, the pandemics of our past, such as SARS in 2003 and flu in 1918-1920, 1957, 1968 and 2009, decayed due to other reasons and not because the viruses evolved to cause milder symptoms 59 .
During the ongoing COVID-19 pandemic it is therefore required to continuously adapt mitigation strategies, drugs, and vaccines in order to mitigate the impact during an early phase of the establishment of this virus. The world has employed massive sequencing programs (based on infected patients, animal, or large-scale sewage monitoring samples) that continuously provide information on changes in the viral genome. While this ensures an early detection of altered genomes, it by principle lacks information about the changed characteristics of the emerging variants. Therefore, only further and detailed investigations like structural analysis and the prediction of the biological impact of occurring mutations allow us to be one step ahead of the virus by helping to forecast altered virulence, transmissibility, or immune evasion potential.
Given the enormous number of emerging genomes, a fast tool such as Catalophore Halos can help to rapidly identify changes in RBD/hACE2 interface fields and guide a decision of which variants to simulate in full atomistic detail. Difference Halos, in particular, provide qualitative information about the changes of specific physico-chemical properties when comparing surfaces of SARS-CoV-2 VOC RBDs. Changes in electrostatics or hydrophobicity, as observed for the Omicron variant (Fig. 2), may be associated with significant alterations in biomolecular interactions and the viral pathogenicity. This can affect viral transmissibility, the binding affinity between RBD and hACE2, immune escape or the range of hosts. Omicron difference Halos indicate a significant shift to positive electrostatic charges at the Omicron RBD surface compared to the other VOCs, which is in line with previous reports 60,61 . Specifically, T478K, Q493R, Q498R, N440K and E484A cause positively charged patches or decrease negative charges at the RBD surface, which is clearly visible in the difference Halos. R493 might interact with negative charges around E35, whereas R498 might be attracted by D38 of hACE2. Further, Table 1. SARS-CoV-2 VOCs relative binding free energies and ratios of dissociation constants in relation to the wild type predicted by our empirical scoring function compared to results determined through surface plasmon resonance by Barton et al. 32 for a subset of these VOCs. www.nature.com/scientificreports/ K417N, occurring in Omicron and Beta RBD, was shown to disrupt the salt bridge to D30 of hACE2 and to neutralizing antibodies, which has also been described for K417T of Gamma. Thus, these mutations could facilitate the immune evasion of these virus variants but weaken the binding to hACE2 50,62 . N501Y was described to be a key residue in hACE2 binding as it favors a hydrophobic cation-π interaction to K353 of hACE2 and a π-π stacking interaction with hACE2 Y41 62 . The local increase of hydrophobicity caused by N501Y is clearly visible in the Alpha, Beta, Gamma and Omicron difference Halos (difference to the wildtype, Fig. 2a). We further examined the difference Halo outcome by focussing on the effects on overall RBD-hACE2 binding affinity. Based on the LIE method, we have developed an empirical binding affinity estimator that predicts pretty accurate binding energies for spike-RBD-hACE2 complexes at moderate wall-clock run times, especially using massive cloud computing facilities. With the necessary caution, this technique helps to raise flags to indicate potential higher infectiousness. The precision of our SARS-CoV-2 spike-RBD-hACE2 binding affinity model is mainly due to the high number of replicates used to achieve a satisfactory convergence in binding energies and both a solid and sufficiently large training set from one source, thanks to Zahradnik et al. 34 With a mean absolute error of less than 2 kJ/mol, our predictions are in the same error range as experimental binding energy methods 63 . To the best of our knowledge, this is the first application of an LIE-based model to the prediction of the binding free energy of an entire protein-protein interface and it is interesting to observe a very uncommon relation between VDW and electrostatic weights. The extremely low binding energy contribution of the latter, already obtained in other published fitting approaches 64 , is obviously caused by the enormous background noise in the case of the huge interface consisting of roughly 500 atoms and a distinct sensitivity of the Coulomb potential to positional displacements causing large fluctuations as published in the supplementary information of a recent article 13 . This uncertainty cannot be compensated by using 50 trajectory replicates and becomes even worse with increasing interface cutoffs or by handling the spike RBD as the ligand. The latter effect is most likely due to the exclusive presence of loops rather than defined secondary structure elements and the absence of hACE2 as a structure-stabilizing binding partner. Fortunately, force field based VDW interaction energies of the "small" 5 Å interface region reveal remarkable stability over the 50 replicates and, on average, correlate well with experimental binding affinities. A reasonable explanation for this stability and correlation is the mathematical short-range nature of the Lennard-Jones potential and, thus, a strong dependence of the binding affinity on the proteins' steric fit. Nevertheless, here as well, the model accuracy strongly suffers under background noise if the interface is specified significantly larger.

Experimental affinities Predicted
The binding affinity is usually quantified by the equilibrium dissociation constant (K D ), which is used to evaluate and classify the strength of biomolecular interactions. The smaller the K D value, the greater the binding affinity of the ligand for its target. The actual K D value relevant to a concrete biological situation depends on the physiological environment, e.g., salt concentration, temperature, or pH. Therefore, both the measured and modeled absolute K D values are only valid within the observation range, making the comparison between unrelated data sources complex. Thus, we reported K D ratios relative to the wild type which clearly show a structural, biologically explainable increase in the binding affinity of all VOCs, particularly pronounced for Alpha and the new VOC Omicron.
Most of our computational results associated with model training as well as application to SARS-CoV-2 VOCs are in agreement with experimental findings. For instance, the relative order of binding affinity predictions for spike VOCs are exactly in line with experimental observations 32,48,50 of WT, Alpha, Beta, and Gamma. Moreover, the outstanding binding strength associated with Alpha (N501Y) as well as the decrease in binding affinity due to the immune escape mutation (K417N/K417T) in addition to N501Y in case of Beta and Gamma is also perfectly reflected. According to our findings, Omicron achieved a substantially higher increase in binding affinity than all other VOCs except for Alpha. Compared to Beta and Gamma, Omicron bears 13 other mutations, six of which were described to enhance binding to hACE2 13,[35][36][37][38]60 . Other mutations within the Omicron RBD were proposed to weaken the binding affinity to hACE2 and further enhance immune escape due to weaker binding to antibodies 60 . This could explain the second strongest binding affinity to hACE2 for Omicron next to Alpha, despite the drastic changes in electrostatics and hydrophobicity observed in the Omicron difference Halos (Fig. 2a). Since the initial submission of this publication in early December 2021, additional experimental binding affinities for Omicron have been published by different groups. In Table S5 we took the opportunity to compare some of these experimental Omicron RBD-ACE2 binding data to our predictions. In contrast to Han et al., who described the Omicron binding affinity to be lower compared to the wild type 60 , other recent publications showed an RBD binding affinity to ACE2 approximately two times higher compared to the wild type [65][66][67] and K D values (VOC RBD to hACE2) from Cameroni et al. 67 correlated well with our predictions. Our initial conclusion still remains true, this new variant requires our special attention and monitoring.
In summary, our ESF derived from the LIE model allows us to estimate high-level binding affinities between the spike RBD and hACE2 regions. However, although we found the model fairly accurate when compared to experimental values, predicted affinities should be interpreted as binding "trends" instead of absolute K D values. As noted elsewhere [16][17][18] , increased binding affinity often leads to increased infectivity of SARS-CoV-2. Given the experience with the global wave of VOC Delta, it is therefore now necessary to closely monitor the spread and impact of Omicron in the coming months. Initial medical reports 68 from South Africa 69 , Botswana, and Europe 70 with a relatively small, non-representative sample of patients revealed a changed clinical picture with a comparatively mild course of the disease 71 . As population immunity increases, either through infection or vaccination, steady modulation of immune-evading mutations might contribute to a permanent establishment of SARS-CoV-2, wherein Omicron apparently shows the potential of playing a significant role.

Methodology
Analysis of spike RBD diversity. For the analysis of RBD diversity, all spike protein sequences that were available by December 6th 2021 at GISAID 21 were downloaded in FASTA format. Processing and analysis of the sequences was performed employing in-house tools in Python. The spike RBD sequences (residues 319-541) were extracted by splitting each spike protein sequence after the amino acid motif "TSNF" and before "NFNG", which refer to the RBD flanking regions. Only sequences consisting of 223 residues (the length of the RBD) were accepted, i.e., insertions and deletions were not considered. Each residue of the retrieved RBD sequences was compared with the respective residue in the reference RBD 21 . Every mismatch except for low-quality or sequencing error residues (indicated by "X") was counted as one mutation. For the analysis of spike RBM diversity, residues 437-508 were considered.
Preparation of spike RBD-hACE2 structures. A 2.47-Å-resolution crystal structure of the SARS-CoV-2 spike RBD bound to hACE2 (PDB-Code: 6m0J) was used as a starting structure. Employing the molecular modeling and MD simulation package Yasara 72 , conformational stress was removed by steepest descent energy minimization followed by simulated annealing (timestep 2 fs, atom velocities scaled down by 0.9 every 10th step). For this purpose, we selected the AMBER14 force field 73 applying an 8 Å force cutoff. SARS-CoV-2 spike-RBD variant sequences were constructed based on the SARS-CoV-2 spike RBD starting sequence from the wild type lineage. Mutant structures were built using homology modeling by implementing Yasara with a maximum of five alignment variations per template and not more than 50 conformations tried per loop. The minimized starting structure served as a template. The final input files contained residues 333-526 of the respective SARS-CoV-2 RBD and residues 19-615 of hACE2 coordinating one zinc ion.

Molecular dynamics simulations.
Again using Yasara, each spike RBD-hACE2 complex was centered in a cuboid simulation box under periodic boundary conditions with a solute-wall distance of 5 Å on every side. This box was filled with explicit solvent molecules of the TIP3P model 74 , approximating 0.997 g/mol water density and with 0.9% sodium chloride ions as well as additional ions for system neutralization. Protein structure topologies and energetics were parameterized according to the Amber14 force field. The system pressure was set to 1 bar using the Manometer1D setting of Yasara, while temperature coupling to 310 K was achieved by velocity rescaling as described by Krieger et al. 75 Long-range electrostatic interactions were treated through Particle Mesh Ewald 76 summation (PME) with an 8 Å cutoff. Our protocol consisted of a steepest descent energy and successive simulated annealing minimization step (following Yasara's standard minimization protocol) and a final production run of 200 ps (respectively 2 ns in case of model optimization) per replicate using a step size of 1 fs for intramolecular and 2 fs for intermolecular forces. Afterwards we extracted intermolecular interaction energy contributions, E vdw and E elec , caused by van der Waals and electrostatic forces, respectively. In addition, hACE2 in its dissociated state underwent the same procedure from simulation box generation up to production MD in order to have both bound and unbound states available for the development of a predictive model.

Development of a binding affinity model.
From the various MD based approaches to binding affinity estimation, the empirical linear interaction energy (LIE) method developed in the 90s by Åqvist et al. 41,42 has provided a remarkable trade-off between predictive accuracy and computational effort 44,63,77,78 . In contrast to complex methods mimicking thermodynamic reaction paths represented by high degrees of decomposition and numbers of long trajectories, it exclusively relies on short simulations of a bound/associated (PL) and an unbound/dissociated state (L) of the ligand as depicted by the embedded graphic in the upper left corner of Fig. 3a. Starting from the linear response approximation 79 , the developers pointed out a strong relationship between the Gibbs free energies of binding on the one side and, on the other side, average (denoted by angle brackets ⟨⟩) differences ΔE vdw and ΔE elec of the ligand's interaction energies with its surrounding atoms, namely protein and solvent atoms in the bound and solely solvent atoms in the unbound case.
The coefficient alpha was calibrated with respect to a small training set of protein-ligand complexes. Subsequent LIE studies revealed significantly better correlations when both coefficients, α and β, were treated as empirical parameters, possibly extended by further features, and fitted to available training sets with known binding affinities [43][44][45] Our final ESF of the spike-hACE2 binding indeed makes use of two empirical parameter weights w vdw and w elec (used instead of α and β and having omitted the angle brackets for the sake of clarity): Since we are dealing with two interacting proteins rather than a protein-ligand system, we tested both proteins independently as ligand molecules simulated in complex as well as freely in solvent. It must also be noted that due to the huge size of these "ligands" and in order to reduce the impact of noise, only interactions of ligand amino acids with at least one atom within a 5 Å (and 10 Å for comparison) environment of the binding partner (e.g. brown colored area of hACE2 object in embedded graphic of Fig. 3a) were evaluated upon energy computation in the bound as well as unbound case.
(1) G = w vdw E vdw + w elec E elec . www.nature.com/scientificreports/ As discussed in the Results and Discussion sections, treating hACE2 as a ligand significantly increased the accuracy of our model, although mutational changes are associated with RBD, whereas the hACE2 sequence remains constant within the training set. Moreover, our model (Eq. 1) treats both energy contributions as empirical parameters fitted to a training set where, in particular, the optimal electrostatic weight deviates enormously from its original value close to 0.5. Since our affinity model deviates from the LIE protocol in several aspects, we decided to call it an empirical scoring function rather than an LIE model, which is nevertheless closely related to it.
For model development and evaluation we employed a sufficiently large training set of 43 spike RBD variants including the wild type (PDB ID 6M0J) in complex with hACE2 for which Zahradnik et al. had recently published K D values obtained through yeast surface display titration 34 . The number of mutations per variant ranges from one to seven, with a maximum distance to hACE2 amounting to 2.5 nm. Using a linear regression model along with least-squares fitting interaction energy weights of the empirical model were fitted to experimental ΔG exp values derived from K D values at a temperature of T = 310 K and using the gas constant R: Gibbs free energies of binding for new variants were then predicted by summing up the two weighted interaction energy contributions according to Eq. (1) and translated back to K D values by inverting Eq. (2).
Our predictive model was validated and optimized through leave-one-out cross-validation, ensuring that none of the 43 variants was inside the training set used to predict its binding energy. Mean absolute errors and squared Pearson correlation coefficients served as measures for model accuracy. In order to achieve a satisfactory binding energy convergence and increase model accuracy, we produced and analyzed 50 replicates of 2 ns MD trajectories and averaged the two interaction energy terms.
Catalophore Halo analysis. Spike RBD Halos were calculated as fields of physico-chemical properties, represented by 3D point clouds covering the entire outer side of the molecular surface of each RBD variant in complex with hACE2. Molecular surfaces were hereby defined with a probe radius of 1.4 Å around the RBD atoms' VDW radii, and the point clouds on the RBD surface had a thickness of 5 Å. These full size Halos were cropped down to the binding interface region, that is the space which is occupied by atoms of the binding partner hACE2 plus corresponding VDW radius plus an additional 5 Å extension around these atoms. It is important to note that hACE2 is not influencing the properties, hence solely the shape, of the spike RBD Halos. The grid spacing of the point clouds was set to 0.75 Å.
Halo point clouds of all spike RBD variants were aligned in space using the iterative closest point algorithm and by taking into account both coordinates and property values. Afterwards, difference Halo point clouds were generated for each pair of two RBD variant Halos under investigation by placing a new point at the geometric center of each pair of two corresponding points of the two primary Halo point clouds. For the definition of pairs of corresponding points we used a distance cutoff of 0.5625 Å. Finally, property values for the new difference cloud points were calculated as the difference of the values associated with the corresponding pair of two primary points, thereby subtracting values of the points in the first primary point cloud from the values of the second primary point cloud. The concept of Catalophore Halos is described in more detail in the supplementary information.
Computational details. The simulation performance was optimized in cooperation with Amazon Web Services (AWS), who supplied the necessary cloud infrastructure in the framework of the diagnostic development initiative. We used clusters of AWS Elastic Computing (EC2) × 86 instances of the C5 und C6i families running 64-bit Amazon Linux 2 with AMI Kernel 5.10, e.g. c6i.8xlarge (32 virtual CPUs), c5.4xlarge (16 vCPUs) and c6i.16xlarge (64 vCPUs) equipped with Intel Xeon Scalable processors. With this configuration the simulation of each variant was completed within 3 to 6 days. Initial proof-of-principle simulations were executed with GROMACS 2021.2 on ARM-based AWS Graviton instances c6g.8xlarge (32 vCPUs).

Data availability
Publicly available datasets were analyzed in this study. This data can be found here: https:// www. gisaid. org/. Input and final structure files as well as Pandas Dataframes of interaction energies exported as Python Pickle files generated within this work are available for download at https:// doi. org/ 10. 6084/ m9. figsh are. 17129 771.